#####################################################################################
#                                                                                   #
#    Scripts for breeding songbird point count data in New Jersey Pine Barrens      #
#                                                                                   #
#                           Created by: Philip M. Coppola                           #
#                                                                                   #
#####################################################################################

##### OPEN POPULATION DISTANCE SAMPLING IN JAGS - see Amundsen et al (2014); Kery & Royle (2016) Ch 9.3; YAMAURA & ROYLE (2017)

# Install Packages
#install.packages(c("R2jags", "mcmcplots"))

# Load Packages
packagelist <- c("R2jags", "mcmcplots")
lapply(packagelist, require, character.only = TRUE)

# Set working directory - update as needed
setwd("E:/My Drive/Phil Coppola Dissertation Research files/All Data/Coppola NJPB Songbird/Outputs/")

bdat <- read.csv("../Combined/birds.csv")

# make a cross-table specified by species, site, visit, and distance
y <- table(bdat$species,bdat$site,bdat$visit,bdat$distance)
y <- y[apply(y,1,sum)>0,,,] 
attributes(y)$dimnames <- NULL

# number of detections against distance classes
apply(y,3,sum)

# habitat information
hab <- read.csv("../Combined/treat.csv")
hab <- hab[,3]

# species information
S <- dim(y)[1]                    # number of species
m <- 30	                          # number of augmented species

# number of sites
nsite <- length(hab)
# number of visits
K <- 3

# Create the distance class data
nD <- 2              # Number of distance classes
B <- 100             # upper limit of counting (maximum count distance)
delta <- B/nD        # bin size or width
mdpts <- seq(delta/2,B,delta) # midpoint distance of bins up to max distance

# arrangement of the augmented data
ay <- array(0,dim=c((S+m),nsite,3,2))  # y of augmented
ay[1:S,,,] <- y

# total detections per site and occasion for individual species
nobs <- apply(ay,c(1:3),sum)

# area used in rarefaction
area <- rep(NA,18)
# patch area required for prediction
for(j in 1:18){
  area[j] <- 0.01+0.2*(j-1)
}
area <- c(c(0.01,0.02,0.04,0.08,0.16,0.32,0.64),area[5:18])
narea <- length(area)
trt <- hab

cat("
model{
  
  # hyperparameters
  
  # treatment effects for single functional group
  for(i in 1:2){	                                                     # number of treatment (1 = control, 2 = treatment; dummy)
    mu.t[i] ~ dnorm(0,0.0001)
    sigma.t[i] ~ dunif(0,5)
    tau.t[i] <- 1/(sigma.t[i]*sigma.t[i])
  }
  
  # detectability (distance sampling parameter)
  mu.s ~ dnorm(0,0.0001)T(1,6)                                               # changed for convergence
  sigma.s ~ dunif(0,5)                                                       # changed for convergence
  tau.s <- 1/(sigma.s*sigma.s)
  
  # emigration rate
  mu.q ~ dnorm(0,0.0001)
  sigma.q ~ dunif(0,5)
  tau.q <- 1/(sigma.q*sigma.q)

  psi ~ dunif(0,1)                                                           # inclusion rate that generates wi

  for(i in 1:(S+m)){

    # detection probability is assumed to be constant among the sites
    l.sigma[i] ~ dnorm(mu.s,tau.s)T(1,6)                                     # log-scale sigma; using truncation; former T(3,5)
    log(sigma[i]) <- l.sigma[i]     	                                       # distance function parameter
    
    # emigration rate
    q[i] ~ dnorm(mu.q,tau.q)
    logit(phi[i]) <- q[i]
    
    # indicator variable whether each species is exposed to sampling or not
    w[i] ~ dbern(psi)				

    for(h in 1:2){	                                                         # treatment effects
      trt.eff[i,h] ~ dnorm(mu.t[h],tau.t[h])
      log(lambda[i,h]) <- trt.eff[i,h]
      dens[i,h] <- lambda[i,h]*phi[i]*w[i] 
    }
    
    # Distance sampling detection probability model
      for(b in 1:nD){
        log(g[i,b]) <- -mdpts[b]*mdpts[b]/(2*sigma[i]*sigma[i])             # Half-normal
        f[i,b] <- ( 2*mdpts[b]*delta )/(B*B)                                # Radial density function
        pi.pd[i,b] <- g[i,b]*f[i,b]                                         # Product Pr(detect)*Pr(distribution)
        pi.pd.c[i,b] <- pi.pd[i,b]/sum(pi.pd[i,1:nD])                       # Conditional probabilities
      }
    
    for(s in 1:nsite){
      for(k in 1:K){
        pdet[i,s,k] <- sum(pi.pd[i,1:nD]) # distance class prob
        pmarg[i,s,k] <- pdet[i,s,k]*phi[i]    # marginal prob      
        
        # model part 3: distance class frequency
        ay[i,s,k,1:nD] ~ dmulti(pi.pd.c[i,1:nD],nobs[i,s,k])
        
        # model part 2: total number of detection
        nobs[i,s,k] ~ dbin(pmarg[i,s,k],M[i,s])
      } # end of k loop
      
      # model part 1: abundance model
      M[i,s] ~ dpois(lambda[i,trt[s]]*w[i])
    } # end of s loop
  } # end of i loop
  
  # derived parameters
  R <- sum(w[1:(S+m)])  # estimating regional species pool size
  
  # rarefaction
  for(i in 1:(S+m)){
    for(j in 1:2){
      for(k in 1:narea){
        DENS[i,j,k] <- (dens[i,j]/3.14159)*area[k]   # calculate per area density first, and multiply the area after (plot area = 3.14159 ha)
        OCCU[i,j,k] <- 1 - exp(-DENS[i,j,k])   # occupancy probability
        UNOCCU[i,j,k] <- 1 - OCCU[i,j,k]       # complement of occupancy probability
      }
    }
    
    for(k in 1:narea){
      TOTALOCCU[i,k] <- 1 - prod(UNOCCU[i,,k]) # probability that the species occurs at least one site
    }
  }
  
  for(i in 1:2){
    for(j in 1:narea){
      rareN[j,i] <- sum(DENS[,i,j])
      rareR[j,i] <- sum(OCCU[,i,j])
    }
  }
  
  for(i in 1:narea){
    meanR[i] <- mean(rareR[i,])
    gammaR[i] <- sum(TOTALOCCU[,i])
  }
} # end of the model
",file="hds.txt")

# initial values
inits <- function(){
  list(M=apply(nobs,c(1:2),sum)+5,
       mu.t=rnorm(2),sigma.t=runif(2),
       mu.s=runif(n=1,min=2,max=6),sigma.s=runif(1),
       mu.q=rnorm(1),sigma.q=runif(1),
       trt.eff=matrix(rnorm((S+m)*2),ncol=2),
       q=rnorm(S+m),
       r=rnorm(S+m),
       l.sigma=runif(S+m)+4,
       w=rep(1,(S+m)),psi=runif(1))
}
params <- c("mu.t","sigma.t",
            "mu.q",
            "sigma.q",
            "phi",
            "mu.s","sigma.s",
            "sigma",
            "trt.eff",
            "p",
            # "lambda",
            "dens",
            "R","psi",
            "rareN","rareR",
            "meanR","gammaR"
)

# MCMC
mi <- 220000
ni <- 11000
nb <-  1000
nt <- 1
nc <- 3

# data
data <- list("ay","nobs","K","nsite","trt","nD","B","mdpts","S","m","area","delta","narea")

library(jagsUI)
out.cam <- autojags(data,inits,params,"hds.txt",n.thin=nt,n.chains=nc,n.burnin=nb,iter.increment=ni,max.iter=mi,parallel=TRUE)
write.table(out$summary,"estimateHDS.csv",sep=",")
